DISCONTINUITY INDUCED BIFURCATIONS OF NON-HYPERBOLIC CYCLES 

IN NONSMOOTH SYSTEMS 

ALESSANDRO COLOMBO^ AND FABIO DERCOLEt 

Abstract. We analyse three codimension-two bifurcations occurring in nonsmooth systems, when a non- 
hyperbolic cycle (fold, flip, and Neimark-Sacker cases, both in continuous- and discrete-time) interacts with one 
of the discontinuity boundaries characterising the system's dynamics. Rather than aiming at a complete unfolding of 
the three cases, which would require specific assumptions on both the class of nonsmooth system and the geometry 
of the involved boundary, we concentrate on the geometric features that are common to all scenarios. We show that, 
at a generic intersection between the smooth and discontinuity induced bifurcation curves, a third curve generically 
emanates tangentially to the former. This is the discontinuity induced bifurcation curve of the secondary invariant 
set (the other cycle, the double-period cycle, or the torus, respectively) involved in the smooth bifurcation. The result 
can be explained intuitively, but its validity is proved here rigorously under very general conditions. Three examples 
from different fields of science and engineering are also reported. 
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1. Introduction. This article deals with the analysis of three particular codimension- 
two bifurcations in nonsmooth systems. Broadly speaking, nonsmooth systems are continuous- 
or discrete-time dynamical systems featuring some kind of discontinuity in the right-hand side 
of their governing equations whenever the system's state reaches a discontinuity boundary. 
More specifically, nonsmooth systems include several classes, e.g., piecewise smooth [9; 11], 
impacting [2], and hybrid [1; 17] systems, which have been largely used in the last decades 
as models in various fields of science and engineering (see references above and therein). 

While methods of numerical continuation allow to easily detect and trace bifurcation 
curves in two-parameter planes, understanding the geometry of bifurcation curves around 
codimension-two points is a key to the construction of complex bifurcation diagrams. In the 
domain of smooth dynamical systems, the unfolding of the most common codimension-two 
points is well known (see, e.g., [15]), and this knowledge is exploited in continuation software 
for the automatic switching among bifurcation branches at these points (see, e.g., [8; 19]). 
The same cannot be said for nonsmooth systems, where, though efficient numerical tools for 
bifurcation analysis are finally starting to appear [6; 23], results are still mostly limited to 
codimension-one cases. A reason for this shortcoming can be found in the fact that nons- 
mooth systems exhibit, along with the standard bifurcations of smooth systems, a great num- 
ber of completely new bifurcations, called discontinuity induced bifurcations, that involve the 
interaction of the system's invariant sets with the discontinuity boundaries. Since the char- 
acteristics of these bifurcations depend critically on both the class of nonsmooth system and 
the geometry of the involved boundaries, the number of possible scenarios is huge and, at the 
moment, truly general results are scarce. It goes without saying that codimension-two cases 
involving simultaneous smooth and discontinuity induced bifurcations, named "type II" in 
[14], are even more numerous, and less understood. 

In this article we analyse type II bifurcations of periodic orbits (limit cycles), that is, 
bifurcations involving a periodic orbit (from now on called the bifurcating cycle) that col- 
lides with a discontinuity boundary while being at the same time non-hyperbolic. Rather than 
aiming at a complete unfolding with reference to a particular class of nonsmooth systems, 
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we concentrate on finding those geometric features that are common to all classes: this is 
accomplished by abstracting our analysis from the nature of the involved boundary. As a 
consequence, our results are incomplete, because they only focus on the geometry of bifurca- 
tion curves around the codimension-two point; on the other hand, they apply more in general 
— a feature that should be welcome in a field where peculiarity seems to be the rule. 

In particular, we show that three codimension-one bifurcation curves generically emanate 
from a type II point in a two-parameter plane. One is the smooth bifurcation curve (fold, flip, 
or Neimark-Sacker), while the other two are the discontinuity induced bifurcations of the bi- 
furcating cycle and of the secondary invariant set involved in the smooth bifurcation (the other 
cycle, the double-period cycle, or the torus, respectively). Then we show that, depending on 
the bifurcation, one or both of these curves are tangent to the smooth bifurcation curve. In- 
deed, in the flip and Neimark-Sacker cases, the bifurcating cycle departs from the image of the 
nonhyperbolic cycle, left frozen in state space, at a linear rate with respect to the bifurcation 
parameter, whereas the distance between the period-two cycle or the torus from such an image 
goes as the square root of the parameter perturbation from the bifurcation. As a consequence, 
locally to the codimension-two point, the perturbation required by the secondary invariant 
set to collide with the discontinuity boundary is quadratic with respect to that required by the 
bifurcating cycle. Similarly, in the fold case, the rate at which both cycles approach the image 
of the nonhyperbolic cycle is proportional to the square root of the parameter perturbation, 
so that the discontinuity induced bifurcation curves are both quadratically tangent to the fold 
curve. These rather intuitive results have been observed in many examples, and proved for 
some specific classes of discontinuous systems (e.g., in [5; 14; 20; 21; 22; 24; 26]). The aim 
of this paper is to provide formal support to the above geometric arguments and to prove their 
validity once and for all under very general conditions. 

The ensuing exposition is set into the framework of grazing bifurcations in continuous- 
time, where the discontinuity boundary is smooth, locally to the point of contact with the 
bifurcating cycle, and the contact occurs tangentially. This allows us to keep the terminology 
as coherent as possible, especially in the lack of a uniform terminology across all classes of 
nonsmooth systems. Nonetheless, the reader will realise that our exposition is general and 
applies to any discontinuity induced bifurcation involving a non-hyperbolic cycle in contin- 
uous time or a non-hyperbolic fixed point in discrete time. In fact, our analysis is based on 
the reduction of the nonsmooth flow to a map which is defined and smooth on one side of 
a boundary, while we do not describe the behaviour of the map on the other side. The rest 
of the analysis is based on the obtained map, as if the problem was originally set in discrete 
time. Thus, in practise, we do not make any assumption on the class of nonsmooth systems 
and on the geometry of the discontinuity boundary. 

We begin by stating the problem, introducing the basic notation, and outlining the steps 
that we follow in the main proofs (Sect. 2); then we proceed with the detailed analysis of the 
three generic grazing bifurcations of non-hyperbolic cycles: the grazing-fold, the grazing- 
flip, and the grazing-Neimark-Sacker (Sects. 3-5 and Appendices). Once casted in discrete 
time, grazing bifurcations are more appropriately called border collisions, and this is the 
name we use in this part of the paper. Then we presents three specific applications (Sect. 6) 
and conclude with some future directions. 

2. The framework of analysis. We consider a nonsmooth autonomous flow x(t) = 
$(x(0), t, a) e R n+1 depending on parameters a e R 2 . Namely, the right-hand side of the 
system's ODEs 



(2.1) x{t)= ^(x(t),r,a) 
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Fig. 2.1: A generic (hyperbolic) limit cycle 7 of the nonsmooth flow $. For some a near 
a = 0, the cycle passes close to, but does not touch, the discontinuity boundary T>, so that 
the resulting Poincare map on V is defined, locally to z, only on one side of the discontinuity 
boundary TL. The boundary H divides V into two regions, respectively composed of points z 
from which the orbit of $ does and does not touch V. 



(here and in the following variables and parameters as subscripts denote differentiation) is 
generically smooth, but characterised by zero- or higher-order discontinuities across some 
discontinuity boundaries Dj, defined as the zero set of suitable smooth functions Di(x, a). In 
particular, we can distinguish three types of discontinuity boundaries (see Fig. 2.1): boundaries 
across which the right-hand side of (2.1) is nonsmooth but continuous, so that orbits always 
cross the boundary (V c in the figure); boundaries across which the right-hand side of (2.1) 
is discontinuous, so that sliding motions are possible (T>s); boundaries where the right-hand 
side of (2.1) is formally characterised by impulsive components, which define an instanta- 
neous state transition (or jump) whenever orbits reach the boundary (T>i). 

Forward solutions of system (2. 1) are composed of smooth segments, each corresponding 
to a smooth orbit terminating at a discontinuity boundary, or to a sliding motion. Smooth 
segments are directly connected at crossing and sliding boundaries, while they are connected 
through state jumps at impacting boundaries. Let 7 be a periodic orbit of system (2.1). In 
Fig. 2.1, 7 is composed of four segments, three smooth (solid) orbits and one sliding motion 
(thick orbit), and is characterised by a single state jump (thick dashed connection). 

Suppose that, when a = 0, the cycle 7 grazes (touches tangentially) a discontinuity 
boundary T>, and no other degeneracies occur on T>c, T>j, and T>g. At the same time, suppose 
that 7 is non-hyperbolic at a = (more precisely, the multipliers are not defined at a = 0, 
but the smooth bifurcation curve is a path to a = on which one real or two complex 
conjugate simple multipliers lie on the unit circle). Introduce a Poincare section V along one 
of the segments of 7, say, e.g., the segment touching T> so that the flow reaches T> after V for 
a = 0. Also introduce a coordinate z e R™ on V such that the intersection z of 7 with V lies 
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at z = for a = 0. Then, locally to (z, a) = (0, 0), the flow $ induces a Poincare map 
(2.2) ZHf(z,a) 

(note that the map may not be invertible, e.g., in the presence of sliding motions). Since we 
do not discuss the type of boundary T>, we limit the definition of F to the values of (z, a) in a 
neighbourhood of (0, 0) for which the orbit originating at z does not touch T>. This introduces 
an (n — 1) -dimensional discontinuity boundary H on the Poincare section V such that F is 
defined and smooth on one side of H. In particular, let 

V = {x: D(x, a) = 0}, H = {z : H(z, a) = 0}, 

and assume, without loss of generality, that the flow $ touches T> tangentially while locally 
remaining on the side D(x, a) < 0, and that F(z, a) is defined for H(z, a) < 0. Then, 
the function H can be constructed as follows (see again Fig. 2.1). Define the n-dimensional 
smooth manifold T of the points where the flow is tangent to the level sets of function D: 

T={x:T(x,a) := {& t (x, 0, a), D x (x, a)) = 0} 

(vector D x (x, a) e R" +1 is orthogonal to the level sets of I? at (x, a) and (•, •) is the standard 
scalar product in R" +1 ). As shown in Fig. 2.1, the (n — l)-dimensional intersection between 
T> and T is transformed, backward in time by the flow, into the discontinuity boundary H. 
Thus, H(z, a) can be defined as the value D(x, a) at the point x at which the flow first 
reaches T (forward in time) from the initial condition corresponding to z on V . 

We can now abandon the continuous-time framework, and focus on map (2.2). For some 
a in a neighbourhood of a = 0, the map is characterised by a fixed point z, with H(z, a) < 
and, for a = 0, the fixed point is non-hyperbolic and lies at the origin z = and on the 
discontinuity boundary H. We investigate the bifurcation curves rooted at a = in the 
parameter plane (a\, a 2 ), by considering separately the three generic cases, namely (I) fold 
(one simple eigenvalue equal to 1, Sect. 3), (II) flip (one simple eigenvalue equal to — 1, 
Sect. 4), and (III) Neimark-Sacker (two simple complex conjugate eigenvalues on the unit 
circle, Sect. 5). 

In each case, we proceed as follows. Locally to (z, a) — (0, 0), we consider the re- 
striction of map (2.2) to a parameter-dependent centre manifold Z c . Let u E R" c represent 
coordinates on Z c , n c = 1 in the fold and flip cases, n c — 2 in the Neimark-Sacker case, 
with u = u(z, a) for each z G Z c and a in a neighbourhood of (z, a) = (0, 0), u(0, 0) = 0, 
and let z = z(u, a) denote the inverse transformation. Restricted to the centre manifold, map 

(2.2) reads 

(2.3) u i ► f(u, a) := u(F(z(u, a), a), a). 

and the discontinuity boundary H is given by the zero-set of the function 

(2.4) h(u, a) := H(z(u, a), a). 

We assume that the three following conditions hold: 

(i) Map (2.3) satisfies, at a — 0, all genericity conditions of the corresponding smooth 
bifurcation (see, e.g., [15]). 

(ii) At a = 0, the centre manifold Z c transversely intersects the discontinuity bound- 
ary H at z = (by continuity the transversality persists near (z, a) = (0, 0), see 
Fig. 2.1). Under this condition, the dynamics of map (2.2) near (z,a) = (0,0) is 
captured by that on the centre manifold. In the coordinate u along the centre mani- 
fold the condition becomes h u (0, 0) ^ 0. 
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(iii) Changing a along the smooth bifurcation curve, the non-hyperbolic fixed point 
crosses the discontinuity boundary transversely. This condition ensures that the 
smooth bifurcation curve intersects the border collision curves in a generic way. 
As a first step, we reduce map (2.3) to a normal form (NF) (the fold, flip, and Neimark-Sacker 
normal forms) through a locally invertible change of variable and parameter, say, v = v(u, a), 
(3 = (3(a), where v(0,0) = 0, (3(0) = 0, and u = u(v,(3), a = a((3) denote the inverse 
transformation. Then, second step, we find the expression of the discontinuity boundary (2.4) 
in the new variables and parameters, i.e., 

(2.5) {v : h NF (v, (3) := h(u(v, (3), a((3)) = 0}. 

Finally, third step, we analyse the interaction of the normal form map 

v i ► f NF (v, (3) := v(f(u(v, (3),a((3)), a((3)) 

with the discontinuity boundary (2.5), and we find local asymptotics for the bifurcation curves 
emanating from a = in terms of (ai, a 2 )- ex P ans i° ns - 

The details of the normal form reduction are reported in appendices A.l, B.l, and C.l, 
while the technicalities on step two are reported in Appendices A. 2, B.2, and C.2. The spe- 
cific analytical form taken by condition (iii) in the fold, flip, and NS cases is respectively 
derived in Appendices A. 3, B.3, and C.3 in terms of both the original coordinates z and in 
the coordinates u in the centre manifold. Finally, some details on step three for the Neimark- 
Sacker case are relegated to Appendix C.4. For simplicity of notation, in the following the 
superscript stands for evaluation at (u, a) — (0, 0) or (v, (3) — (0, 0). 

3. Case I: Border-fold bifurcation. Let the dynamics in the centre manifold Z c be 
described by the one-dimensional system 

(3.1) u f(u,a), u e R 1 , 

with f° = (fixed point condition) and /° = 1 (fold condition). Under condition (i), map 

(3.1) can be reduced to normal form (first step, see Appendix A.l) with invertible changes of 
variable and parameter v = v(u, a), (3 — (3(a), becoming 

(3.2) v ^ (3i + v + sv 2 + 0(v 3 ), 

where s — sign(/° u ). In these variables, the fold curve has equation (3\ — in the plane 
((3\- l (3'i), and the corresponding non-hyperbolic fixed point is located at v = 0. 

We now turn our attention to the discontinuity boundary (2.5) (second step, see Appendix 
A.2). Condition (ii), ensuring transversal intersection of the centre manifold Z c and the 
discontinuity boundary H, implies local existence and uniqueness of a smooth function 

a((3)=*l(3 1+ al(3 2 + 0(\\(3\\ 2 ), 

such that the intersection of H with Z c is located at v — cr((3). Then by condition (iii) (see 
Appendix A. 3 for the analytical expression) we know that moving along the fold curve, that 
is, along the /3 2 -axis, the fixed point at v = crosses H at (3 2 = 0. As a consequence, we 
have <r° 2 £ 0. 

We are now ready to find the equation of the border collisions in the plane ((3\,fli) 
(third step). The two fixed points of the normal form map (3.2) are located at w ± (/3) = 
±^/~s(3\ + 0(\\(3\\ 2 ) (v~ being stable and v + unstable for s = 1, and viceversa for s = —1), 
and lie on the discontinuity boundary (2.5) along the curves 



(3.3) 



± = 4a + 4a + o(\\/3\\ 2 ). 
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Fig. 3.1: Border-fold bifurcation. Bifurcation curves: LP, fold (limit point, red); BC S , border 
collision of the stable fixed point (y~ , left; v + , right) of map (3.2) (green), BC", border 
collision of the unstable fixed point (v + , left; v~, right) of map (3.2) (blue). Region labels: 
0, no fixed point in V~(j3) := {v : h NF (v, (3) < 0}; 1, v~ is the only fixed point in V~((3)\ 
2, both fixed points vr 1 lie in V~ ((3). 



Since cr2 7^ 0, equation (3.3) for small \\(3\\ becomes 
(3.4) ±yf^fc~o O p a 2 , 

and gives the asymptotics, locally to j3 = 0, of the two border-collision bifurcation curves 
involving the fixed points w ± . The invertible parameter change (3 = (3(a) easily provides the 
asymptotics in the original a parameters. 

Depending upon the sign of s in the normal form map (3.2), of er2 in (3.4), and of 
in (ii), there are eight generic cases, two of which are reported in Fig. 3.1. The other six can 
be reduced to these two by suitable parameter changes. In fact, the four cases with a^ 2 < 
are symmetric with respect to the /3i-axis to the corresponding cases with <t° 2 > 0, while the 
four cases with ft," < can be reduced to cases with > by changing the sign of s and 
rotating the figure. Note that only half of the /^-axis can be said to belong to the fold curve 
(LP), since along the other half the two fixed points w ± collide at v = on the undescribed 
side of the discontinuity boundary (2.5), i.e., ft NF (0, (3) > 0. 

4. Case II: Border-flip bifurcation. Let the dynamics in the centre manifold Z c be 
described by the one-dimensional system 

(4.1) uw/(a,a), iteR 1 , 

with f° — (fixed point condition) and /£ = — 1 (flip condition). Through a parameter- 
dependent translation, we can ensure that /(0, a) = 0, i.e., that u — is a fixed point for all 
a in a neighbourhood of a = 0. Under condition (i), map (4. 1) can be reduced to normal form 
(first step, see Appendix B.l) with invertible changes of variable and parameter v — v(u, a), 
[3 = (3(a), becoming 



(4.2) 



v I-*- -(1 + 0i)v + sv 3 + 0(v A ), 
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with s = sign((l/4)(/° u ) 2 + (l/6)/° uu ). In these variables, the flip curve has equation 
j3\ = in the plane (fii, f3 2 ), and the corresponding non-hyperbolic fixed point is located at 
v = 0. Moreover, parameters can be chosen so that the border collision of the fixed point in 
the origin has equation /3 2 = 0. 

We now turn our attention to the discontinuity boundary (2.5) (second step, see Appendix 
B.2). Condition (ii), ensuring transversal intersection of the centre manifold Z c and the 
discontinuity boundary H, implies local existence and uniqueness of a smooth function 

a(P)=alp 1 +alp 2 + Om\ 2 ), 

such that the intersection of H with Z c is located at v = <r{(3). Moreover, thanks to the 
parameter choice in (4.2), cr^ = 0, since the fixed point v — lies on H when /3 2 = 0. Then 
by condition (iii) (see Appendix B. 3 for the analytical expression) we know that moving along 
the flip curve, that is, along the /3 2 -axis, the fixed point at v = crosses H at /3 2 = 0. As a 
consequence, we have a^ 2 ^ 0. 

We are now ready to find the equation of the border collisions in the plane (/3i , /3 2 ) (third 
step). Near (v, (3\) — (0, 0) the normal form map (4.2) iterated twice has one fixed point in 
v = (which is also a fixed point of map (4.2)) and two others in w ± (/3) = ±\/s/3i+0(||/3|| 2 ) 
(period-two cycle). In particular, lie on discontinuity boundary (2.5) along the curves 

(43) ±^7p~ 1 =a%l3 2 +Om\ 2 )- 
Since er° 2 ^ 0, equation (4.3) for small \\[3\\ becomes 

(4.4) iv^fr^jfcft. 

and gives the asymptotics, locally to (3 = 0, of the border-collision bifurcation curves involv- 
ing the two points u ± of the period-two cycle. The invertible parameter change (3 — (3(a) 
provides the asymptotics in the original a parameters. 

Depending upon the sign of s in the normal form map (4.2), of ct° 2 in (4.4), and of in 
(ii), there are eight generic cases. However, again, only two cases are relevant (see Fig. 4.1), 
because all others can be reduced to these two by suitable parameter changes. Here, both the 
four cases with cr° 2 < and those with h° a < 0, are symmetric with respect to the /3i-axis to 
the corresponding cases with (T° 2 > or > 0. Also note that only half of the /3 2 -axis can 
be said to belong to the flip curve (PD), since along the other half the fixed point v = lies on 
the undescribed side of the discontinuity boundary (2.5), i.e., ft, NF (0, 0) > 0. Similarly, only 
one of the two branches in (4.4) constitutes the border-collision curve involving the period- 
two cycle (stable, BC|; unstable, BC 2 ), since along the other branch ft, NF (w ± , (3) > 0. 

5. Case III: Border-Neimark-Sacker bifurcation. Let the dynamics in the centre man- 
ifold Z c be described by the two-dimensional system 

(5.1) u^/(ii,a), u e R 2 , 

with f° = (fixed point condition) and with eigenvalues A and A (the overbar stands for 
complex conjugation) of the 2 x 2 Jacobian /° given by 

\(a) = (l + g(a))e if> W, 



with g° = (Neimark-Sacker, NS, condition). As in the flip case, assume that /(0, a) = 
for all a in a neighbourhood of a = 0. Under condition (i), map (5.1) can be reduced to 
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Fig. 4.1: Border-flip bifurcation. Bifurcation curves: PD, flip (period doubling, red); BC^' U , 
border collision of the fixed point v = (stable and unstable branches, blue); BC^", border 
collision of the stable or unstable period-two cycle. Region labels: 0, no fixed point or period- 
two cycle in V~((3) := {v : /i NF (i> , 0) < 0}; 1, v = Ois a fixed point in V~((3) and there is 
no period-two cycle, or it does not lie entirely in V~ (/?); 2, the fixed point v = coexists in 
V~ (f3) with the period-two cycle. 



normal form in polar coordinates (first step, see Appendix C.l) with invertible changes of 
variable and parameter p = p(u, a),ip — ip(u, a), (3 — (3(a), becoming 



(5.2a) 
(5.2b) 



p h+ p(l + ft + a(/V) + P 4 R(P, <P, (3), 
^^ + 0(a(/3)) +/9 2 Q(p,^/?), 



where a ^ 0. In these variables, the NS curve has equation ft — in the plane (ft , ft), 
and the corresponding non-hyperbolic fixed point is located at v — (with v\ = Ke(pe l>p ) 
and t>2 = Im(pe i¥J )). Moreover, parameters can be chosen so that the border collision of the 
fixed point in the origin has equation ft = 0. 

We now turn our attention to the discontinuity boundary (2.5) (second step, see Appendix 
C.2). Condition (ii), ensuring transversal intersection of the centre manifold Z c and the 
discontinuity boundary H, implies local existence and uniqueness of a smooth function 

a{(3) =^ 1 ft+^ 2 ft + 0(||/5|| 2 ), 



measuring the distance between the origin and the boundary, with positive/negative values 
if /i NF (0,ft) is negative/positive, in order to make u(f3) differentiable at (3 = 0. Moreover, 
thanks to the parameter choice in (5.2), er^ = 0, since the fixed point v = lies on TC when 
ft — 0. Then by condition (iii) (see Appendix C.3 for the analytical expression) we know 
that moving along the NS curve, that is, along the ft -axis, the fixed point at v = crosses H 
transversely at ft = 0. As a consequence, we have cr^ 2 ^ 0. 

We are now ready to find the equation of the border collisions in the plane (ft , ft) (third 
step). Near (3 — 0, the normal form map (5.2) has an fixed point in p — and a closed 
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Fig. 5.1: Border-NS bifurcation. Bifurcation curves: NS, Neimark-Sacker (red); BC SU , bor- 
der collision of the fixed point v = (stable and unstable branches, blue); GR S ", grazing 
of the stable or unstable torus (green). Region labels: 0, no fixed point or invariant curve in 
V~(0) := {v : /i NF (i>, (3) < 0}; 1, v = Ois a fixed point in V~((3) and there is no invariant 
curve, or it does not lie entirely in V~(/3); 2, both the fixed point v = and the invariant 
curve lie in V~{(3). 
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invariant curve that is contained in the annular region 

{«"•"> : f^) {l -^ l,2) s " 5 + «" /2) ' * e |0 ' 2 " ] } • I < ■> < 1 

(5.3) 

(see Appendix C.4). The two circles delimiting the annular region (5.3) touch the discontinu- 
ity boundary along the curves 



(5-4) (1 ± pj- 1 ' 2 ) = + 0(\\(3\n 

Since cr^ 2 ^ 0, equation (5.4) for small ||/3|| becomes 

(5-5) y^J^a° 2 /3 2 , 

and gives a unique asymptotic, locally to j3 = 0, for the grazing bifurcation curves of both 
circles. The same asymptotic therefore holds for the grazing bifurcation involving the invari- 
ant curve (the uniqueness of the bifurcation curve is granted by the elliptical shape of the 
invariant curve near (3 = 0). Again, the invertible parameter change (3 — (3(a) provides the 
asymptotic s in the original a parameters. 

Depending upon the sign of a in the normal form map (5.2) and of a°a in (5.5), there 
are four generic cases. However, again, only two cases are relevant (see Fig. 5.1), because 
those with <r° 2 < are symmetric with respect to the /3i-axis to the cases with <r° 2 > 0. 
Also note that only half of the /?2-axis can be said to belong to the NS curve, since along the 



10 



A. Colombo and F. Dercole 



other half the fixed point v = lies on the undescribed side of the discontinuity boundary 
(2.5), i.e., /i NF (0, (3) > 0. Similarly, only half of the parabola in (5.5) constitutes the grazing 
bifurcation curve involving the invariant curve (stable, GR S ; unstable, GR U ), since along the 
other half the invariant curve is composed of points v with h NF (v, [}) > 0. 

6. Examples. We now present three specific examples, one for each of the three codimension- 
two bifurcations analysed in the previous sections. The three examples deal with different 
classes of nonsmooth systems (an impacting, a hybrid, and a piecewise smooth system) and 
describe interesting applications in different fields of science and engineering (ecology, social 
sciences, and mechanics). 

An impacting model of forest fires. For an example of border-fold bifurcation, we 
consider the forest fire impacting model presented in [7; 18]. The model describes the vege- 
tational growth with the following two (smooth) ODEs: 



one for the surface layer (bush, B) and one for the upper layer (trees, T). Fire episodes are 
represented by instantaneous events (impacts), that occur when the biomasses (B, T) of the 
two layers reach one of three specified impacting boundaries: a bush ignition threshold p B K B 
triggering bush-only fires that map the bush biomass to X B p B K B , < X b ,Pb < 1; a tree 
ignition threshold p T K T triggering trees-only fires that map the trees biomass to X T p T K T , 
< \t,Pt < 1; and the segment connecting points (ctbKb, PtKt) an d (pbK b , &tKt), 
< <jb < Pb, < ut < Pt, triggering mixed fires with post-fire conditions suitably 
assigned as a function of pre-fire conditions (see [18] for more details). 

For the parameter setting n = 0.375, = 0.0625, a = 0.43, K B — Kt = 1, 
Pb = 0.85, pt = 0.93, Xb = 0.03, At = 0.01, o B = 0.61, ot = 0.3 (corresponding 
to Mediterranean forests), the system is characterised by a globally stable period-one cycle 
composed of a growth orbit and a mixed fire. Numerical continuation (by means of AUTO07P 
[10]) of the cycle in the parameter plane (p B , pt) identifies two (codimension-one) bifurca- 
tions: a fold (red curve in Fig. 6.1) and a grazing of the growth orbit with the bush ignition 
threshold (blue curve). The two curves merge together at the border-fold bifurcation (black) 
point and, as predicted by the analysis carried out in Sect. 3, the grazing bifurcation of the 
unstable cycle involved in the fold (green curve) emanates tangentially to the fold curve from 
the codimension-two bifurcation point. 

A hybrid model of two-party democracies. For an example of border-flip bifurca- 
tion, we consider the hybrid model presented in [3] for describing the dynamics of two-party 
democracies. The model describes the evolution of the size of two lobbies (of sizes Lp and 
Lr), one associated to each party (parties D and R, respectively), and assumes that the in- 
dividuals belonging to the lobby of the party at the government erode the welfare (W) at a 
rate proportional to the size of the lobby; a lobby can grow only as long as its party is at the 
government, and decays otherwise; a small fraction of the lobbyists not at the government 
defect and switch to the other lobby; elections are held once every T years, and people vote 
for the party that has the less damaging lobby at the time of the elections. Altogether, the 
dynamics is captured by two sets of ODEs, namely 




W = r(l- W-a D L D )W 1 
Ld = (e D a D W - d D )L D + k R L R , 
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Lr = (— djj - k R )L R , 

when the £>-party is at the government, and 

W = r(l - W-a R L R )W, 
L D = {—do - k D )L D , 
L R = (e R a R W - d R )L R + k D L D , 

when the i?,-party is at the government. Here, r is the intrinsic growth rate of the welfare, a 
represents the aggressiveness of a lobby, e is the recruitment coefficient of a lobby, and d and 
k are respectively the rate at which individuals abandon the lobbies or defect. In the region 
of the state space where ao^D < clrL r (cldLd > clrL r ) the D-lobby (i?-lobby) is less 
damaging and thus wins the elections. The condition anLu = a R L R therefore defines the 
discontinuity boundary (see [3] for more details). 

In the (an, T) plane, with parameters a R = 1, r = 0.2, en = e R = 6, do = d R = 1.8, 
kD = k R = 0.06, the system has a very complex bifurcation diagram (see for example 
Fig. 1 in [3]). In particular, near cid = 0.38, T = 3.2, a flip (red curve in Fig. 6.2) and a 
border collision (blue curve) of a period-2T cycle meet at the border-flip (black) point and, 
as predicted by the analysis carried out in Sect. 4, a border collision of the period-4T cycle 
(green curve) emanates from the codimension-two point tangentially to the flip curve. 

A piecewise smooth model of railway wheelset dynamics. For an example of border- 
NS bifurcation, we consider a two degrees of freedom piecewise smooth model of a sus- 
pended railway wheelset with dry friction dampers, subject to a sinusoidal disturbance repre- 
senting the deformations of the track. The model is based on that presented in [13; 25], where 
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the track deformation was not taken into account, and its analysis will be published elsewhere. 
Since a detailed explanation of the equations and parameters goes beyond the scope of this 
paper, here we only report the equations and describe a few key parameters (see [13] and [25] 
for the details). The model consists of the following piecewise smooth equations: 

&\ = x 2 , 

±2 = — {~2F X - 2K s x 1 - sign(x 2 )At), 

TO 

i 3 = x 4 , 

± 4 = ~(-2AF y ), 

where 

X\ = Xx + asin(a->t), x 2 = x 2 + aa>cos(u;i), 
fi = — sech(ai 2 )) + /Lt s sech(ai 2 )), 

p & F r p tyFr _ J £ r C (l - ^ + ^) if C£ r < 3 Mt , 

* (fit otherwise, 

Here w = 27r V/ Z, a and / are the amplitude and wavelength of the sinusoidal disturbance, V 
is the speed of the wheelset, and A measures the conicity of the wheels. The system's state 
space is therefore partitioned in four regions, depending on the signs of x 2 and of C£ r — 3/Zt, 
so that x 2 = and C£ r = 3/ij define two discontinuity boundaries. 
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The system's dynamics was studied, with TC-HAT [23], in the (V, A) plane, with the 
following values of the parameters: m = 1022, K s = le6, 1 = 678, A = 0.75, a = 0.001, 
fj, d = 1000, a = 50, fx s = 1200, * = 0.54219, $ = 0.60252, C = 6.5630e6, fi t = le5, 
tq = 0.4572, I = 10. For large values of V, a grazing of a stable cycle with the boundary 
X2 = and a NS take place (blue and red in Fig. 6.3), and meet at the border-NS (black) 
point. Then, by systematically evaluating 1000 iterations (after transient) of the Poincare 
map of the torus on a suitable cross-section, and by continuing the line on which the obtained 
torus image grazes the discontinuity boundary induced on the cross-section, we were able to 
trace an approximation of the grazing curve of the torus (green in Fig. 6.3). More rigorous 
methods, based for example on discretisation of the invariant curve (see, e.g., [4; 12]) could 
be used to obtain a more precise estimate of the quadratic coefficient. This lies however 
beyond of the scope of this paper. As predicted by the analysis carried out in Sect. 5, the 
curve emanates from the codimension-two point tangentially to the NS curve. 

7. Concluding remarks. We have analysed the geometry of bifurcation curves around 
three codimension-two bifurcations in nonsmooth systems, namely the border-fold, the border- 
flip, and the border-Neimark-Sacker. Rather than aiming at the complete unfolding of the 
dynamics of a particular class of nonsmooth systems (e.g., piecewise smooth, impacting, 
or hybrid) dealing with a particular geometry of the involved discontinuity boundary (e.g., 
smooth or corner), we have focused on those results which are general to all scenarios. Our 
approach applies to continuous-time as well as discrete-time systems, and basically consists 
of the analysis of a discrete-time (Poincare) map defined only on one side of a boundary in its 
state space. Explicit genericity conditions are listed and explained for each codimension-two 
case. 

Of course, the weakness of this approach is that it cannot provide the complete unfolding 
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of the bifurcation, but its power resides in its generality: as shown in the three examples that 
we have reported, it applies to a very broad class of nonsmooth systems and it may be relevant 
in various fields of science and engineering. 

The natural sequel of this work would certainly aim at more detailed results, and possibly 
at the complete unfolding, of the codimension-two bifurcations analysed here, with specific 
reference to some smaller class of nonsmooth systems. 
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Appendix A. Border-fold bifurcation. In the case of the border-fold bifurcation, con- 
ditions (i-iii) in Sec. 2, expressed in the variable u of the centre manifold, are summarised 
below: 

(i-a) fl u + 0, 

(i-b) /° ± 0, 

(ii) K + 0, 

(iii) f° h" f - h" f° f° ^ f° h° f° - h° f° f° 

V"V J uu a±Ja 2 uJ ua\J a 2 ~ J uu' v a 2 J ai 'u J ua 2 J ai 

Note that (i.b) is redundant, since it is implied by (iii). 

A.l. Step one. To reduce map (3.1) to normal form we follow [15], where however 
a e R, while here a e R 2 . The variable change v ~ v(u,a) is formally the same as in 
[15], while parameter change that we use is (3 = /3(a) = \a(^i(a))\^i(a), n\(a) = f® a a\ + 

^ 2 («) = -./5 , ct2 « 1 +./5' c , 1 « 2 +c , (ll«ll 2 ), - / 2 («(m))+o(|| q; (^)||), 

with a(0) = (l/2)/° u ^ because of (i.a). The inverse transformations have the following 
derivatives: 



9 f° 9 

7/° - ?;° - —A rv° S° - ^ a - 

U v — I rO I ' U (il — °a a f3 2 i °a ~ M > a 2 ~ | r ill y^O 

I J uu I J uu I J uu II I •/ a 



_f0 

J a 2 

f° 

J Qi 



A.2. Step two. Consider the discontinuity boundary (2.5). The variable and parameter 
change v — v(u,a), (3 = (3(a) is invertible near (u,a) = (0,0), so that condition (ii) 
implies that ft,^ F (0,0) = h^u® ^ 0, i.e., local existence and uniqueness, by the implicit 
function theorem, of a smooth function 

a(f3)=alP 1+ alf3 2 + 0(\\pf), 

such that h NF (a/3, (3) = for small 1 1 (3\ | , so that the intersection of the discontinuity boundary 
H with the centre manifold Z c is located at v — a((3). 

We now prove, using condition (iii), that <r^ 2 ^ 0. By differentiating both sides of 
h NF (a(P), (3) — 0, i.e., of h(u(a((3), (3), a([3)) = 0, with respect to /3 2 , taking into account 
the derivatives in Appendix A.l, and evaluating at (3 = we get 



fr°^° 2 + h l a% _L / / ,0 _ h lfu ai \ ,0 _ I h _ h u.fua 2 \ f0 

02 h° u U v h u \\f a f {{ ^ f uu J^ 2 fL ) U J 

Thanks to (i)-(iii), this ensures that <jp 2 ^ 0. 

A.3. Genericity conditions (ii) and (iii). In the original coordinates z of map (2.2), 
condition (ii) requires H^t>° ^ 0, where v is the unit eigenvector of F z associated to the 
eigenvalue 1. 

Consider now the fold curve defined by the system 

F(z,a)-z = 0, 
(A.l) F z (z,a)v-v = 0, 

(v, v)-\ = 0. 

In the space (z, v, a), condition (iii) means that the tangent vector to the fold curve is not 
tangent to the surface 

(A.2) H(z,a) = 
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at (z, a) — (0, 0). The tangent vector to the fold curve is the null vector of the Jacobian of 
(A.l), so that bordering such Jacobian with the linearisation of (A. 2) and imposing that the 
resulting square matrix is nonsingular at (z, v, a) — (0, v°, 0), i.e., 



det 



( FS-I 





ai 


F° 


\ 




F°-I 


F° v° 


F° u° 





2(V>) T 















H° 


H° 

a 2 


) 



7^0, 



we impose that the fold curve (A.l) intersects the surface (A.2) transversely, i.e., condition 
(iii). This is nothing but requiring that the system (A.l), (A.2) be regular at (z, v, a) = 
•:<>./'".•>:•. 

Equation (A.l), restricted to the centre manifold, becomes 



f(u,a)-u = 
f u (u,a)-l = 

and by the same reasoning, we obtain the condition 



dct 



'0 i 

u 


fO 

J Oil 


f° 


f° 


fO 

J uai 


fO 


K 


h° 


h° 



which is equivalent to (iii) since /" = 1 (fold condition). 

Appendix B. Border-flip bifurcation. In the case of the border-flip bifurcation con- 
ditions (i-iii) in Sec. 2, expressed in the variable u of the centre manifold, are summarised 
below: 

(i-a) \{fL) 2 + \fuuu^^ 
(i.b) fl a ? 0, 
(ii) K ? 0, 

v 111 / Jtioi"a 2 I Jua 2 "ai 

Note that (i.b) is redundant, since it is implied by (iii). 

B.l. Step one. Once again, to reduce map (4.1) to normal form, we use the same vari- 
able change v = v(u, a) as in [15], while the parameter change is j3\ = j3\(a) = g(ai, a.2), 
02 = 02(a) = h(0,a), with f u (0,a) = —(1 + g{a)). The inverse transformations have 
derivatives 



1 



1 



withc° = (l/4)(/0j 2 + (l/6)/l 



f° h° - f° h° 

Jua\' b a.2 J UOL2 OL\ 

^ because of (i.a). 



-f° 

J u 



f° 



B.2. Step two. Consider the discontinuity boundary (2.5). As in the border-fold case, 
the variable and parameter change v — v(u, a), (3 — [3(a) is invertible near (u, a) — (0, 0), 
so that condition (ii) implies that h^ F (0, 0) 7^ and, by the implicit function theorem, that 
the intersection of the discontinuity boundary H with the centre manifold Z c is located at 



v = a(0)=a o fjl /3 1 +a o j2 + O 



for some smooth function a. 
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The parameter change obviously makes cr^ = 0. We now prove that ct° 2 ^ 0. By 
differentiating both sides of /i NF (ct(/3) , /3) = 0, i.e., of h(u(a(/3), /?), a(/3)) = 0, with respect 
to /?2, taking into account the derivatives in Appendix B.l and evaluating at j3 2 = we get 



h° u u% 



+ h° a a% 



h?.u°„ 



h° 



where condition (iii) ensures that hP a ap 2 — 1. Thus (i)-(iii) imply that crj^ 7^ 0. 

B.3. Genericity conditions (ii) and (iii). In the original coordinates z of map (2.2), 
condition (ii) requires H®v° ^ 0, where v is the unit eigenvector of F z associated to the 
eigenvalue —1. 

Consider now the flip curve defined by the system 



(B.l) 



F(z,a)-z = 0, 
F z (z,a)u + u = 0, 
{v, v)-\ = 0. 



Similarly to the border-fold case, condition (iii) is equivalent to 



dot 










F° 

ai 


F° 


F° + I 


F° u° 


F° v 




x zai 


zai 


2v T 











H° 


H° 

02 



^0. 



Equation (B. 1), restricted to the centre manifold, becomes 

/(it, a) — u = 0, 
f u (u,a) + l = 0. 

Proceeding along the same lines we obtain the condition 

1 f° f° 

J Ctl J Oil 




" 

UOL\ 



/; 
hi 



f° 

J ua.2 

h° 



^0, 



which is equivalent to (iii) since /° 
condition). 



(/(0, a) = by assumption) and /° = -1 (flip 



Appendix C. Border-Neimark-Sacker bifurcation. In the case of the border-NS bi- 
furcation, conditions (i-iii) in Sec. 2, expressed in the variables u of the centre manifold, are 
summarised below: 

(i.a) e tk9 ° 1 for k = 1, 2, 3, 4, 

(i.b) the first Lyapunov coefficient of the NS normal form (a , see later) is nonzero, 
(i-c) 9°a + 0, " 

(ii) K ± 0, 

(iii) g° ai h° a ^ 9° a2 h° ai 

Note that (i.c) is redundant, since it is implied by (iii). 

C.l. Step one. Once again, to reduce map (5.1) to normal form, we use the same vari- 
able change w = w(u,a) (with w = v\ + W2) as in [15], while the parameter change 
(3 = /3(a) is formally the same as in B.l. The inverse transformations u = u(w,w,(3) and 
a = a(f3) have derivatives 



u w {0,0,0) = q°, u*b(0,0,0) 



5°, 



^(0,0,0) = 0, aj 



ya 1 n a 2 Ba 2 ' l a 1 



i)a\ 
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Fig. C.l: A. Local representation of the discontinuity boundary E (thick line) for small \\v\\ 
and 1 1/3 1 1 as a straight (dashed) line tangent to E in the point of minimum distance of E from 
the origin v = (case with <r(f3) > 0). Since \\(3\\ is small, the direction ip m of minimum 
distance is close to the direction iph of vector ft,^ F (0, 0). For <p <G (ipo, ipi) (shaded area), 
the discontinuity boundary E can be represented in coordinates (r, ip). B. The annular region 
(5.3) (shaded area) containing the invariant curve (thick closed line) of the normal form map 
(5.2) and the (dashed) circle approached by the invariant curve as ft — > 0. 



C.2. Step two. Denote by E the discontinuity boundary (2.5), where v E R 2 . Again, 
the variable and parameter change v — v(u,a), (3 — /3(a) that we used is invertible near 
(u,a) — (0,0), so that condition (ii) implies that h^ F (0, 0) = hP u v? v ^ 0, where now 
h^ F (0, 0) and ft," are in R 2 (row vectors) and w° is a 2 x 2 nonsingular matrix. Geometrically, 
see Fig. C.1A, this means that for small \\v\\ and \\(3\\ we can represent the discontinuity 
boundary (2.5) as a straight line almost orthogonal to h^ F (0, 0) and slightly displaced from 
v = in the direction of /i^ F (0, 0). 

Let ifih be the angle of vector h^ F (0, 0) with respect to axis v\. Technically, 

VH = arctan 27r (/j* F (0, 0), /i* 2 F (0, 0)), 

where arctan 27r is the four-quadrant inverse tangent in [0, 2ir]. For any in a neighbour- 
hood of ifih, introduce axis r passing from the origin v = with direction ip, so that positive 
and negative values of r measure the distance from the origin along directions ip and ip ± tt, 
respectively (see Fig. C. 1 A). Coordinates (r, p) are like polar coordinates, but allow differen- 
tiation with respect to r at r = 0. We can therefore express the discontinuity boundary (2.5) 
as 



where 



E = {(r >¥ >) :^((rcos(^),rsin(^)),/3)=0}, 



d_ 

dr 



h {(rcos(tp h ),rsm(tp h )),0) 



= C(o,o) 



r=0 



cos((p h ) 
sm(ip h ) 



^0 



(recall that, by definition of ph, hy F (0, 0) is proportional to (cos(iph), sin(<^^))), so that, by 
the implicit function theorem, we can represent E explicitly as r — 5(ip, 0), 5(p, 0) = 0, for 
some smooth function S defined for p in an open neighbourhood (ip , pi) of ip^. 
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Now, define <p m ((3) := axgmin ve( ^ 0iVl) {|5(<p,/3)|}for/3 ^ and note that lim^o <Pm (P) = 
<Ph> so that we can set (p® n = Lp h . Then, the minimum distance of £ from the origin v = is 
given by the absolute value of 

a{fi) := % m (/3), 0) = *lf3i + + 0(||/3|| 2 ), 

while its sign says whether the minimum is realised along the direction <p m {(3), if positive, or 
(fim(P) i 7r, if negative. In the first case (see Fig. C.1A), v = is a fixed point of the normal 
form map (5.2), since /i NF (0, (3) < 0, while v — lies on the undescribed side of X in the 
second case, i.e., /i NF (0, /?) > 0. 

Similarly to the border-flip case, the parameter change implies that cr2 = 0. We now 
show that er^ 7^ 0. By differentiating both sides of /i NF ((<5(<p, (3) cos(y>), <5(<p, /?) sin(yj)), /?) = 
0, i.e., of 

&(«(%, /?)e*%,/3)e-*/3),a(/3)) = 0, 

with respect to /3 2 , taking into account the derivatives in Appendix C.l, and evaluating at 

ft = we get 

, , n , hlu 02 {Q,W) + hla% 1 



/#(u°e*+u° e -*) 2^Re( g °e^)' 

which is well defined for p ~ Lp h thanks to (ii). Indeed, vP w e %iph + u%e~ ttph is nothing but 
d/dr(u(re liph ,re~ liph , 0))| r= o and thus gives the direction of u-perturbations from u = 
corresponding to r-perturbations from r — along the direction ipt, so that, by definition of 
iph, Re(q°e liph ) is proportional to h°. Finally, we have 



(recall that 5(ip, 0) — for all <p e (<£o, <£i)), so that op 2 ^ thanks to conditions (ii) and 
(iii) (which is necessary to show that h^a^ = 1). 

Note that, in order to evaluate ct° 2 , we need an expression for <ph in terms of variables u. 
For this we can write u as a function of (v, (3), i.e., 



u = u(v,(3) = u(vi + iv2,vi - iv2,(3) 
(u must be read as a function of (w, w, (3) in the right-most side), so that 

< = u w (0, 0,0) + ^(0,0,0) - 2Rc(g°), 

< = u^O, 0,0)* -«*((), 0,0)i = -2Im(<7°), 

and 

Vh, = arctan 27r (hlu° i7 hlu° 2 ) = arctan 27r (/i°Re(g°), -/i"lm(g )) . 

C.3. Genericity conditions (ii) and (iii). Condition (ii) requires H° (Rc(nu°), Im(^ )) ^ 
0, where is the complex unit eigenvector of F z associated to the eigenvalue (1 + g)e t6 '. 
The NS curve is described by the system 

F(z,a)-z = 0, 
(L j 9(a) = 
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where, for any given a, g(a) e R is obtained by solving the system 

F z (Q,a)v - (1 + g)e ie v = 0, 
(!/,!/> -1 = 0, 
Re(z/) T Im(i/) = 0, 

in the variables (<?, 6, v). In the space (z, a) condition (iii) means that the tangent vector to 
the NS curve is not tangent to the surface 

H(z,a) = 

at (z, a) = (0, 0). Similarly to the border-fold and -flip cases, condition (iii) is equivalent to 




Equation (C.l), restricted to the centre manifold, becomes 

f(u,a)-u = 0, 
g(a) = 0. 

By the same reasoning we obtain the condition 

( fu-I fl 

det g° u 5 ° ff ° I ± 



"0 T 
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f° 


f° 
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Q° 




h° 
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which is equivalent to (iii) since /° = (/(0, a) — by assumption) and /° — I is nonsin- 
gular (condition (i.a) (k = 1)). 

C.4. Step three. In this appendix we show that near (i = the closed invariant curve of 
the normal form map (5.2) is contained in the parameter-dependent annular region (5.3) (we 
adapt the material from [16], Chap. 5). 

Assume the supercritical case, i.e., a < 0, so that the invariant curve exists for (3\ > 
and is stable. The annular region shrinks around the circle of equation 



(C2) p= ^__|_ ^€[0,277], 

with 0(/?7)-width (see Fig. C.1B) and map (5.2a) maps p into p + Ap with Ap = p(/3i + 
a((3)p 2 +p 3 R(p^,P)) and 



> P (2f3-; +1/2 - aV + o{pT)) if o < P < ^^(1 - /?r 1/2 ) ; 

< p(- 2 /t 1/2 - pv + oiP 3 / 2 )) if P > f^a + pr 1/2 ). 



Ap < 



Thus the orbits of map (5.2) enter the annular region if 7 < 1 (the term 0l +1 ^ 2 dominates the 
others and determines the sign of Ap), so that with 1/2 < 7 < 1 the stable invariant curve 
remains in the annular region for small \\/3\\. Similarly, in the subcritical case, a > 0, the 
invariant curve exists for (3\ < and is unstable, and the orbits of map (5.2) exit the annular 
region if 7 < 1. Again, with 1/2 < 7 < 1, the invariant curve remains in the annular region 
for small ||/3||. 



